libraries

library(tidyverse)
library(brms)
library(bayesplot)
library(here)
library(glue)
library(scales)

theme_set(theme_light())

source(glue("{params$common_dir_str}/brms_model.R"))
source(glue("{params$model_dir_str}/model_prior.R"))

data

obs_only <- 
  read_csv(glue("{params$model_dir_str}/data/stimulation_obvs.csv")) %>%
  mutate(subj = as_factor(subj),
         obs_degree = error,
         error = obs_degree * (pi/180))
## Parsed with column specification:
## cols(
##   subj = col_double(),
##   subj_index = col_double(),
##   stimulation = col_double(),
##   error = col_double()
## )

full model

print(bprior_full)
##              prior class        coef group resp   dpar nlpar bound
## 1 normal(3.8, 0.4)     b   intercept            circSD            
## 2   normal(0, 0.4)     b stimulation            circSD            
## 3  normal(0, 0.25)    sd   Intercept  subj      circSD            
## 4  normal(0, 0.25)    sd stimulation  subj      circSD            
## 5   normal(0, 1.5)     b   intercept             theta            
## 6     normal(0, 1)     b stimulation             theta            
## 7   normal(0, 0.5)    sd   Intercept  subj       theta            
## 8   normal(0, 0.5)    sd stimulation  subj       theta
iter = 4000
warmup = 2000
cores = 4
chains = 4
n_post_samples = (iter - warmup) * chains

fit_full <- brm(bform_full, obs_only, prior = bprior_full, 
                family = vm_uniform_mix, stanvars = stanvars,
                warmup = warmup, iter = iter, cores = cores, chains = chains, 
                control = list(adapt_delta = 0.99), inits = 0, 
                file =  glue("{params$save_dir_str}/fit_full"))
## Compiling the C++ model
## Start sampling

fit check

divergences

model_fit <- fit_full

#check neff and rhat and divergences
np <- nuts_params(model_fit)
rhat <- brms::rhat(model_fit)
neff_rat <- neff_ratio(model_fit)

np %>% 
  filter(Parameter == "divergent__") %>%
  summarise(n_div = sum(Value))
##   n_div
## 1     0

rhat

mcmc_rhat(rhat) + yaxis_text(hjust = 1) + scale_x_continuous(breaks = pretty_breaks(6))
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

neff ratio

mcmc_neff(neff_rat) + yaxis_text(hjust = 1)

trace plots

mcmc_trace(as.array(model_fit$fit))

reduced model: no pre to post change on circSD parameter

print(bprior_DcircSD_null)
##              prior class        coef group resp   dpar nlpar bound
## 1 normal(3.8, 0.4)     b   intercept            circSD            
## 2  normal(0, 0.25)    sd   Intercept  subj      circSD            
## 3  normal(0, 0.25)    sd stimulation  subj      circSD            
## 4   normal(0, 1.5)     b   intercept             theta            
## 5     normal(0, 1)     b stimulation             theta            
## 6   normal(0, 0.5)    sd   Intercept  subj       theta            
## 7   normal(0, 0.5)    sd stimulation  subj       theta
fit_DcircSD_null <- brm(bform_DcircSD_null, obs_only, prior = bprior_DcircSD_null, 
                        family = vm_uniform_mix, stanvars = stanvars,
                        warmup = warmup, iter = iter, cores = cores, chains = chains, 
                        control = list(adapt_delta = 0.99), inits = 0, 
                        file =  glue("{params$save_dir_str}/fit_DcircSD_null"))
## Compiling the C++ model
## Start sampling

fit check

divergences

model_fit <- fit_DcircSD_null

#check neff and rhat and divergences
np <- nuts_params(model_fit)
rhat <- brms::rhat(model_fit)
neff_rat <- neff_ratio(model_fit)

np %>% 
  filter(Parameter == "divergent__") %>%
  summarise(n_div = sum(Value))
##   n_div
## 1     0

rhat

mcmc_rhat(rhat) + yaxis_text(hjust = 1) + scale_x_continuous(breaks = pretty_breaks(6))
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

neff ratio

mcmc_neff(neff_rat) + yaxis_text(hjust = 1)

trace plots

mcmc_trace(as.array(model_fit$fit))

reduced model: no pre to post change on pMem parameter

print(bprior_DpMem_null)
##              prior class        coef group resp   dpar nlpar bound
## 1 normal(3.8, 0.4)     b   intercept            circSD            
## 2   normal(0, 0.4)     b stimulation            circSD            
## 3  normal(0, 0.25)    sd   Intercept  subj      circSD            
## 4  normal(0, 0.25)    sd stimulation  subj      circSD            
## 5   normal(0, 1.5)     b   intercept             theta            
## 6   normal(0, 0.5)    sd   Intercept  subj       theta            
## 7   normal(0, 0.5)    sd stimulation  subj       theta
fit_DpMem_null <- brm(bform_DpMem_null, obs_only, prior = bprior_DpMem_null, 
                      family = vm_uniform_mix, stanvars = stanvars,
                      warmup = warmup, iter = iter, cores = cores, chains = chains, 
                      control = list(adapt_delta = 0.99), inits = 0, 
                      file = glue("{params$save_dir_str}/fit_DpMem_null"))
## Compiling the C++ model
## Start sampling

fit check

divergences

model_fit <- fit_DpMem_null

#check neff and rhat and divergences
np <- nuts_params(model_fit)
rhat <- brms::rhat(model_fit)
neff_rat <- neff_ratio(model_fit)

np %>% 
  filter(Parameter == "divergent__") %>%
  summarise(n_div = sum(Value))
##   n_div
## 1     0

rhat

mcmc_rhat(rhat) + yaxis_text(hjust = 1) + scale_x_continuous(breaks = pretty_breaks(6))
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

neff ratio

mcmc_neff(neff_rat) + yaxis_text(hjust = 1)

trace plots

mcmc_trace(as.array(model_fit$fit))

reduced model: no pre to post change on pMem or circSD parameter

bprior_DcircSD_DpMem_null
##              prior class        coef group resp   dpar nlpar bound
## 1 normal(3.8, 0.4)     b   intercept            circSD            
## 2  normal(0, 0.25)    sd   Intercept  subj      circSD            
## 3  normal(0, 0.25)    sd stimulation  subj      circSD            
## 4   normal(0, 1.5)     b   intercept             theta            
## 5   normal(0, 0.5)    sd   Intercept  subj       theta            
## 6   normal(0, 0.5)    sd stimulation  subj       theta
fit_DcircSD_DpMem_null <- 
              brm(bform_DcircSD_DpMem_null, obs_only, prior = bprior_DcircSD_DpMem_null, 
                  family = vm_uniform_mix, stanvars = stanvars,
                  warmup = warmup, iter = iter, cores = cores, chains = chains, 
                  control = list(adapt_delta = 0.99), inits = 0, 
                  file = glue("{params$save_dir_str}/fit_DcircSD_DpMem_null"))
## Compiling the C++ model
## Start sampling

fit check

divergences

model_fit <- fit_DcircSD_DpMem_null

#check neff and rhat and divergences
np <- nuts_params(model_fit)
rhat <- brms::rhat(model_fit)
neff_rat <- neff_ratio(model_fit)

np %>% 
  filter(Parameter == "divergent__") %>%
  summarise(n_div = sum(Value))
##   n_div
## 1     0

rhat

mcmc_rhat(rhat) + yaxis_text(hjust = 1) + scale_x_continuous(breaks = pretty_breaks(6))
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.

neff ratio

mcmc_neff(neff_rat) + yaxis_text(hjust = 1)

trace plots

mcmc_trace(as.array(model_fit$fit))

loo compare

#need to do this only once?
expose_functions(fit_full, vectorize = TRUE)

fit_full <- add_loo(fit_full, file = glue("{params$save_dir_str}/fit_full"))
fit_DcircSD_null <- add_loo(fit_DcircSD_null, file = glue("{params$save_dir_str}/fit_DcircSD_null"))
fit_DpMem_null <- add_loo(fit_DpMem_null, file = glue("{params$save_dir_str}/fit_DpMem_null"))
fit_DcircSD_DpMem_null <- add_loo(fit_DcircSD_DpMem_null, file = glue("{params$save_dir_str}/fit_DcircSD_DpMem_null"))

print("fit_full")
## [1] "fit_full"
fit_full$loo
## 
## Computed from 8000 by 504 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -805.2 15.9
## p_loo         6.7  0.2
## looic      1610.4 31.7
## ------
## Monte Carlo SE of elpd_loo is 0.0.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
print("fit_DcircSD_null")
## [1] "fit_DcircSD_null"
fit_DcircSD_null$loo
## 
## Computed from 8000 by 504 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -805.5 15.8
## p_loo         6.3  0.1
## looic      1610.9 31.6
## ------
## Monte Carlo SE of elpd_loo is 0.0.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
print("fit_DpMem_null")
## [1] "fit_DpMem_null"
fit_DpMem_null$loo
## 
## Computed from 8000 by 504 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -804.5 15.9
## p_loo         6.1  0.1
## looic      1608.9 31.8
## ------
## Monte Carlo SE of elpd_loo is 0.0.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
print("fit_DcircSD_DpMem_null")
## [1] "fit_DcircSD_DpMem_null"
fit_DcircSD_DpMem_null$loo
## 
## Computed from 8000 by 504 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -804.9 15.8
## p_loo         6.0  0.1
## looic      1609.9 31.6
## ------
## Monte Carlo SE of elpd_loo is 0.0.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.
loo_compare(fit_full, fit_DcircSD_null, fit_DpMem_null, fit_DcircSD_DpMem_null, criterion = "loo")
##                        elpd_diff se_diff
## fit_DpMem_null          0.0       0.0   
## fit_DcircSD_DpMem_null -0.5       0.8   
## fit_full               -0.8       0.4   
## fit_DcircSD_null       -1.0       1.0